##R Markdown
library(leaflet)
popup = c("Robin", "Jakub", "Jannes")
leaflet() %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%
addMarkers(lng = c(-3, 23, 11), #経度
lat = c(52, 53, 49), #緯度
popup = popup)
library(leaflet)
popup2 = c("Las Vegas", "London", "Switzerland", "Sydney")
leaflet() %>%
addProviderTiles("OpenStreetMap.HOT") %>%
addCircleMarkers(lng = c(-115, 0, 9, 151), #経度
lat = c(36, 51, 47, -34), #緯度
popup = popup2)
install.packages("sf")
## Installing package into '/usr/local/lib/R/site-library'
## (as 'lib' is unspecified)
## Warning in install.packages("sf"): installation of package 'sf' had non-zero
## exit status
install.packages("raster")
## Installing package into '/usr/local/lib/R/site-library'
## (as 'lib' is unspecified)
install.packages("spData")
## Installing package into '/usr/local/lib/R/site-library'
## (as 'lib' is unspecified)
#remotes::install_github("Nowosad/spDataLarge")
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
library(spData)
library(spDataLarge)
world
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>
names(world)
## [1] "iso_a2" "name_long" "continent" "region_un" "subregion" "type"
## [7] "area_km2" "pop" "lifeExp" "gdpPercap" "geom"
summary(world["lifeExp"])
## lifeExp geom
## Min. :50.62 MULTIPOLYGON :177
## 1st Qu.:64.96 epsg:4326 : 0
## Median :72.87 +proj=long...: 0
## Mean :70.85
## 3rd Qu.:76.78
## Max. :83.59
## NA's :10
plot(world)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all
world["lifeExp"] #これの属するdata.frameはsfの一部なので、必ずgeomが付いてくる
## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 2
## lifeExp geom
## <dbl> <MULTIPOLYGON [°]>
## 1 70.0 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01…
## 2 64.2 (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.…
## 3 NA (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.68729…
## 4 82.0 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.…
## 5 78.8 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05…
## 6 71.6 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048…
## 7 71.0 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 …
## 8 65.2 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.27…
## 9 68.9 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1…
## 10 76.3 (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, …
## # … with 167 more rows
class(world)
## [1] "sf" "tbl_df" "tbl" "data.frame"
plot(world["lifeExp"])
world %>% dplyr::select(lifeExp)
## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 2
## lifeExp geom
## <dbl> <MULTIPOLYGON [°]>
## 1 70.0 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01…
## 2 64.2 (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.…
## 3 NA (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.68729…
## 4 82.0 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.…
## 5 78.8 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05…
## 6 71.6 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048…
## 7 71.0 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 …
## 8 65.2 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.27…
## 9 68.9 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1…
## 10 76.3 (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, …
## # … with 167 more rows
world %>% dplyr::select(lifeExp) %>% st_drop_geometry() #st = spatial type
## # A tibble: 177 x 1
## lifeExp
## * <dbl>
## 1 70.0
## 2 64.2
## 3 NA
## 4 82.0
## 5 78.8
## 6 71.6
## 7 71.0
## 8 65.2
## 9 68.9
## 10 76.3
## # … with 167 more rows
#st_drop_geometry(world["lifeExp"]) 上と同じ文
world_asia = world[world$continent == "Asia", ]
plot(world_asia)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all
asia = st_union(world_asia)
## although coordinates are longitude/latitude, st_union assumes that they are planar
plot(asia)
plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red") #重畳表示。前のやつでreset=FALSEとし、次のでadd=TRUEとする。動かない。
#chapter2-2 is vector
# the rbind function simplifies the creation of matrices
## MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
multi_point = st_multipoint(multipoint_matrix)
multipoint_matrix
## [,1] [,2]
## [1,] 5 2
## [2,] 1 3
## [3,] 3 4
## [4,] 3 2
multi_point
## MULTIPOINT ((5 2), (1 3), (3 4), (3 2))
class(multi_point)
## [1] "XY" "MULTIPOINT" "sfg"
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)
## LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)
lnd_point = st_point(c(0.1, 51.5)) # sfg=simple feature geometry object
lnd_geom = st_sfc(lnd_point, crs = 4326) # sfc object
lnd_attrib = data.frame( # data.frame object
name = "London",
temperature = 25,
date = as.Date("2017-06-21")
)
lnd_attrib
## name temperature date
## 1 London 25 2017-06-21
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom) # sf object, ここでようやく使えるデータになる
lnd_sf
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
## Geodetic CRS: WGS 84
## name temperature date geometry
## 1 London 25 2017-06-21 POINT (0.1 51.5)
#chapter 2-3 is raster
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") #spDataパッケージのデータの読み出し
new_raster = raster(raster_filepath)
new_raster #raster layer class. ラスター1つのクラス
## class : RasterLayer
## dimensions : 457, 465, 212505 (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333 (x, y)
## extent : -113.2396, -112.8521, 37.13208, 37.51292 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /usr/local/lib/R/site-library/spDataLarge/raster/srtm.tif
## names : srtm
## values : 1024, 2892 (min, max)
plot(new_raster)
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
r_brick = brick(multi_raster_file)
r_brick
## class : RasterBrick
## dimensions : 1428, 1128, 1610784, 4 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 301905, 335745, 4111245, 4154085 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs
## source : /usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif
## names : landsat.1, landsat.2, landsat.3, landsat.4
## min values : 7550, 6404, 5678, 5252
## max values : 19071, 22051, 25780, 31961
plot(r_brick) # 1~4は光の波長ごとの成分。青、緑、赤、近赤外
plotRGB(r_brick)
# ndviは植生の豊かさを示す。近赤外と赤から算出できる。
ndvi = (r_brick[[4]] - r_brick[[3]]) / (r_brick[[4]] + r_brick[[3]])
plot(ndvi)
library(sf)
library(raster)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(stringr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
##
## extract
library(spData)
sel_area = world$area_km2 < 10000
summary(sel_area)
## Mode FALSE TRUE
## logical 170 7
small_countries = world[sel_area, ]
small_countries
## Simple feature collection with 7 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -67.24243 ymin: -16.59785 xmax: 167.8449 ymax: 50.12805
## Geodetic CRS: WGS 84
## # A tibble: 7 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 PR Puerto Ri… North Am… Americas Caribbean Depe… 9225. 3534874 79.4
## 2 PS Palestine Asia Asia Western … Disp… 5037. 4294682 73.1
## 3 VU Vanuatu Oceania Oceania Melanesia Sove… 7490. 258850 71.7
## 4 LU Luxembourg Europe Europe Western … Sove… 2417. 556319 82.2
## 5 <NA> Northern … Asia Asia Western … Sove… 3786. NA NA
## 6 CY Cyprus Asia Asia Western … Sove… 6207. 1152309 80.2
## 7 TT Trinidad … North Am… Americas Caribbean Sove… 7738. 1354493 70.4
## # … with 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>
plot(small_countries["pop"])
world %>%
filter(area_km2 < 10000) %>%
dplyr::select(pop) %>%
plot()
##GDPperCapitaがtop10の国をピックアップ
world %>%
top_n(n=10, wt=gdpPercap) %>%
dplyr::select(pop) %>%
plot()
#ワースト10
world %>%
top_n(n=-10, wt=gdpPercap) %>%
dplyr::select(pop) %>%
plot()
#日本はGDPcapitaで上から何番目か
world %>%
top_n(n=22, wt=gdpPercap) %>%
top_n(n=-1, wt=gdpPercap) %>%
dplyr::select(gdpPercap) %>%
plot()
#22番目
#解答例
world %>%
filter(name_long == "Japan") %>%
dplyr::select(gdpPercap) %>%
st_drop_geometry() %>%
as.numeric()
## [1] 37337.32
gdp_sort = world %>%
arrange(desc(gdpPercap)) %>%
dplyr::select(name_long, gdpPercap) %>%
st_drop_geometry()
gdp_sort
## # A tibble: 177 x 2
## name_long gdpPercap
## * <chr> <dbl>
## 1 Qatar 120860.
## 2 Luxembourg 93655.
## 3 Brunei Darussalam 76089.
## 4 Kuwait 70832.
## 5 United Arab Emirates 63943.
## 6 Switzerland 57218.
## 7 United States 51922.
## 8 Saudi Arabia 49958.
## 9 Ireland 48898.
## 10 Netherlands 45668.
## # … with 167 more rows